rm(list=ls())

load("TSCS_data.RData")

# Figure E.8

# linear

var_list <- c("lag(polarization,1)","lag(clientelism_gov_seat_share,1)","lag(region_citizen_support,1)",
              "lag(region_antipluralism_gov_seat_share,1)","lag(region_liberal_democracy,1)","lag(democratic_stock,1)",
              "lag(GDP_capita,1)","lag(life_expectancy,1)")

n <- length(var_list)

id1 <- unlist(lapply(1:n,function(x)combn(1:n,x,simplify=F)),recursive=F)

f1 <-lapply(id1,function(x)
  paste("citizen_support ~ lag(citizen_support,1:2)  + antipluralism_gov_seat_share_lag +",paste(var_list[x],collapse="+")))

models_lm <- lapply(f1,function(x) plm(as.formula(x),
                                       data=panel,
                                       model = "within",
                                       effect = "twoways"))

coeffs_lm <- vector()
se_lm <- vector()

for (i in 1:255) {
  coeffs_lm[i] <- models_lm[[i]][["coefficients"]][["antipluralism_gov_seat_share_lag"]]
}

vcm_lm <- lapply(models_lm, function(x) plm::vcovBK(x, cluster=c("group")))
rse_lm <- lapply(vcm_lm, function(x) sqrt(diag(x)))

for (i in 1:255) {
  se_lm[i] <- rse_lm[[i]][["antipluralism_gov_seat_share_lag"]]
}

coeff_df <- as.data.frame(cbind(coeffs_lm,se_lm))
coeff_df %>% mutate(margin_of_error = se_lm*1.65,
                    CI_lower = coeffs_lm - margin_of_error,
                    CI_upper = coeffs_lm + margin_of_error) %>%
  arrange(coeffs_lm) %>%
  mutate(model = row_number())->
  coeff_df

pdf(file = "fig_e8a.pdf",width = 5,height = 9)

ggplot(coeff_df,aes(x=coeffs_lm,y=model))  + 
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), colour="grey") +
  xlab("Coefficient") + ylab("Model") +
  geom_vline(xintercept = 0) + geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

dev.off()

# GMM

n<-length(var_list)
id1<-unlist(lapply(1:n,function(x)combn(1:n,x,simplify=F)),recursive=F)

f1 <-lapply(id1,function(x)
  paste("citizen_support ~ lag(citizen_support,1:2)  + antipluralism_gov_seat_share_lag +",paste(var_list[x],collapse="+"),"|lag(citizen_support, 3:5)"))

models_gmm <- lapply(f1,function(x) pgmm(as.formula(x),
                                         data=panel,
                                         effect = "individual",
                                         model = "onestep",
                                         transformation = "ld",
                                         indexes = c("country", "year"),
                                         robust = T))

coeffs_gmm <- vector()
se_gmm <- vector()

for (i in 1:255) {
  coeffs_gmm[i] <- models_gmm[[i]][["coefficients"]][["antipluralism_gov_seat_share_lag"]]
}

vcm_gmm <- lapply(models_gmm, function(x) plm::vcovHC(x, cluster=c("group")))
rse_gmm <- lapply(vcm_gmm, function(x) sqrt(diag(x)))

for (i in 1:255) {
  se_gmm[i] <- rse_gmm[[i]][["antipluralism_gov_seat_share_lag"]]
}

coeff_df <- as.data.frame(cbind(coeffs_gmm,se_gmm))
coeff_df %>% mutate(margin_of_error = se_gmm*1.65,
                    CI_lower = coeffs_gmm - margin_of_error,
                    CI_upper = coeffs_gmm + margin_of_error) %>%
  arrange(coeffs_gmm) %>%
  mutate(model = row_number())->
  coeff_df

pdf(file = "fig_e8b.pdf",width = 5,height = 9)

ggplot(coeff_df,aes(x=coeffs_gmm,y=model)) + 
  geom_errorbarh(aes(xmin=CI_lower, xmax=CI_upper), colour="grey") +
  xlab("Coefficient") + ylab("Model") +
  geom_vline(xintercept = 0) + geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

dev.off()
